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Abstract 

Thylacinidae is an extinct family of Australian and New Guinean marsupial carnivores, comprizing 12 known species, the 
oldest of which are late Oligocene (~24 IVla) in age. Except for the recently extinct thylacine {Thylacinus cynocephalus), most 
are known from fragmentary craniodental material only, limiting the scope of biomechanical and ecological studies. 
However, a particularly well-preserved skull of the fossil species Nimbacinus dicksoni, has been recovered from middle 
Miocene (~ 1 6-11 .6 Ma) deposits in the Riversleigh World Heritage Area, northwestern Queensland. Here, we ask whether N. 
dicksoni was more similar to its recently extinct relative or to several large living marsupials in a key aspect of feeding 
ecology, i.e., was N. dicksoni a relatively small or large prey specialist. To address this question we have digitally 
reconstructed its skull and applied three-dimensional Finite Element Analysis to compare its mechanical performance with 
that of three extant marsupial carnivores and T. cynoceplialus. Under loadings adjusted for differences in size that simulated 
forces generated by both jaw closing musculature and struggling prey, we found that stress distributions and magnitudes 
in the skull of N. dicksoni were more similar to those of the living spotted-tailed quoll {Dasyurus maculatus) than to its 
recently extinct relative. Considering the Finite Element Analysis results and dental morphology, we predict that N. dicksoni 
likely occupied a broadly similar ecological niche to that of D. maculatus, and was likely capable of hunting vertebrate prey 
that may have exceeded its own body mass. 
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Introduction 

Thylacinids first appear in the Australian fossil record during 
the late Oligocene (~24 Ma) and include the largest representa- 
tives of the Dasyuromorphia, i.e., families Thylacinidae, Dasyur- 
idae and Myrmecobiidae [1-8]. A wide range of feeding ecologies 
are known within the order. They include omnivores, insectivores, 
small prey specialists, hypercarnivores and osteophageous species 
[9]. Variation in the dentition, skuU shape and body size (~1- 
60 kg) of thylacinids suggests considerable trophic diversity within 
the family [10,1 1]. In addition to the recently extinct thylacine or 
Tasmanian 'tiger' {Thylacinus cynocephalus), eleven extinct species of 
thylacinid have been described [4,12-18]. Up to five species may 
have co-existed in the Riversleigh World Heritage Area, north- 
western Queensland between the late Oligocene (~24 Ma) to 
middle Miocene (16-11.6 Ma) [12,16]. See Table SI for the 
temporal and geographic distribution of all thylacinid species. 

The Riversleigh thylacinids inhabited forests [19,20]. These 
regions were also occupied by an assortment of other carnivorous/ 



omnivorous taxa, including 'giant' carnivorous rat-kangaroos 
[Ekaltadeta spp.), crocodiles (Mekosuchinae, e.g. Bam darrowi and 
Trilophosuchus rackhami), flightless dromornithid birds, marsupial 
lions (Thylacoleonidae), bandicoots (Peramelemorphia), dasyurids 
(Dasyuridae), pythons (Pythonidae), madtsoiid snakes and the 
world's oldest known venomous snakes [21-23]. Subsequent 
drying of the Australian continent from the late Miocene (11.6— 
5.3 Ma) led to the gradual replacement of forest environments 
with open woodlands, shrublands and grasslands [19-21]. These 
changes appear to broadly correlate with declining thylacinid 
diversity [24]. 

To date, interpretations of the ecology and feeding behavior of 
fossil thylacinids have been largely qualitative. This is, at least in 
part, because most extinct species are known only from jaw 
fragments and teeth. The near-complete skull oi Nimbacinus dicksoni 
[4,14,18], a medium-sized thylacinid, provides an opportunity to 
more fully investigate feeding ecology in an extinct thylacinid. 

Mimhacinus dicksoni was approximately 5 kg in body mass [10]. 
Fossils of jV'. dicksoni have been recovered from Oligocene-Miocene 
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(~24— 5.3 Ma) deposits in the Riversleigh World Heritage Area, 
northwestern Queensland and Bullock Creek, Northern Territory 
[12,14,16,18]. Its dentition is less specialized than that of the 
species of Thy/acinus, but broadly similar to that of the living 
dasyurid, the spotted-tailed quoU [Dasyums maculatus) in the 
arrangement and geometry of molar shearing crests typically 
associated with camivory [18,25]. 

To date conflicting evidence has been presented regarding the 
body size of prey N. dicksoni may have hunted. Predictions of bite 
force adjusted for body mass, based on applicati(m of 2D beam 
theor)'-, have suggested that N. dicksoni may ha\x" tak(;n relatively 
large prey, as does the sHghdy smaller D. maculatus [26]. However, 
shape analysis of the cranium has suggested that the species may 
have been restricted to smaller prey and/or included a higher 
proportion of invertebrate food in its diet [1 1]. 

The loads imposed on an animal during prey acquisition and 
feeding play an important role in the evolution of its skull 
morphology [27]. Testing hypotheses regarding the relationship 
between the form and function of skulls from extinct species 
requires an understanding of this relationship in living animals 
[28]. A comparative biomechanics approach involving living 
analogues has increasingly been applied to predict the feeding 
ecology and predatory behavior of extinct species [29-34]. To 
gain further insight in the feeding ecology of jV! dicksoni, here we 
perform a biomechanical analysis of the skull of N. dicksoni to 
predict its mechanical behavior in response to loads simulating the 
capture and processing of prey. 

Finite Element Analysis (FEA) is a computer modeling approach 
now commonly used by biologists and paleontologists to examine 
and compare mechanical performance in biological structures in 
comparative contexts [35-41]. In FEA, continuous structures, 
such as the skull, are divided into discrete, finite numbers of 
elements, allowing the prediction of mechanical behavior for 
complex geometric shapes. The structure is analyzed in the form 
of a matrix algebra problem that is solved with the aid of a 
computer [42]. 

Studies of feeding ecology for thylacinids have primarily focused 
on the most recently extinct member of the family, T. cynocephalus, 
which survived in Tasmania until 1936 [43]. Our understanding 
of the ecology of T. cynocephalus is chiefly based on morphological 
comparisons and 2D beam theory [11,26,44,45], as weU as 
anecdotal accounts of their behavior in the wild [46,47]. Elbow 
joint morphology of T. cynocephalus evidently most closely resembles 
that of extant ambush predators, a compromise between efficient 
distance locomotion and the ability to manipulate and grapple 
with prey [48] . Three-dimensional biomechanical modeling of the 
skuU of T. cynocephalus, extant dasyurids and an introduced 
Australian predator {Canis lupus dingo) have suggested potential 
limitations in prey body size [32,35]. 

Morphological and biomechanical comparisons including sym- 
patric nativ(' predators in Tasmania indicat(; that the diet of T. 
cynocephalus may have overlapped considerably with the two largest 
extant marsupial carnivores, the Tasmanian devil (Sarcophilus 
harrisii) and spotted-tailed quoU {D. maculatus) [35,44,49]. These 
species represented the three largest marsupial carnivores in 
Tasmania at the time of European settlement. 

Sarcophilus harrisii is the largest living marsupial carnivore (mean 
adult male weight 8.7 kg, female 6.1 kg) [50]. They are the only 
specialized scavengers among living marsupials, filling a broadly 
similar ecological niche to that of osteophagous hyenas [44]. They 
are also opportunistic hunters and are known to prey on mammals 
that may exceed their body mass [51,52]. Relative to its body size, 
its predicted bite force is greater than that of any other extant 
mammal studied to date [26]. 



QuoUs are represented by four extant species in Australia and 
two in New Guinea [53]. The largest quoU, Dasyums maculatus 
(maximum weight 7 kg) has a broad diet mainly consisting of 
mammals and insects, but will occasionally feed on birds and 
reptiles [54,55]. Larger prey species constitute a higher proportion 
of the diet of adult male D. maculatus, while females and immature 
D. maculatus more frequendy feed on smaUer bodied mammals and 
invertebrates [56,57]. The northern quoU {Dasyums hallucatus) is the 
smallest and most arboreal of the four Australian quoUs, weighing 
up to 1.2 kg [58,59]. Although primarily insectivorous, this active 
hunter can feed on a variety of foods: fruits, smaU mammals, birds, 
reptiles, frogs and carrion [60-62]. 

In this study we aim to determine whether jV. dicksoni was 
capable of IdlUng large prey relative to their body size, or was 
restricted to catching relatively small bodied species. By digital 
reconstruction of its skuU and applied 3D FEA, we compare its 
mechanical performance with that of three extant marsupial 
carnivores; S. harrisii, D. maculatus and D. hallucatus. We use 
previously applied scaling procedures [31] to account for 
diflerences in body mass, allowing for comparison of results 
between species. We predict that N. ditksoni will show similar 
distributions and magnitudes of craniomandibular stress to that of 
D. maculatus due to similarities in their dental morphology and size 
[18]. We also include T. cynocephalus to establish whether the 
biomechanical performance of N. dicksoni more closely resembles 
that of this larger, more derived thylacinid than dasyurids. A 
general assessment of the phylogenetic relationships of dasyur- 
omophians, including taxa examined in this study is shown in 
Figure S 1 . We test if the biomechanical patterns and the inferred 
feeding aspects in T. cynocephalus in a recent study [35] are derived. 
We hypothesize that the relatively long rostrum of T. cynocephalus 
win result in higher stresses in the skull during biting and prey 
procurement than JV. dicksoni and other dasyuromorphians. 

Materials and Methods 

Specimens 

The skull oi Nimbacinus dicksoni (QMF36357) was recovered from 
AL90 Site on the Gag Plateau of the Riversleigh World Heritage 
area in northwestern Queensland. Precise locality details can be 
provided on application to the Queensland Museum. This 
specimen was collected under permits issued by Queensland 
Department of Environment and Heritage and Environment 
Australia, and registered in the paleontological coUections of the 
Queensland Museum, Brisbane, Australia. All necessary permits 
were obtained for the described study, which comphed with aU 
relevant regulations. 

The biomechanical performance of the Mimbacinus dicksoni skuU 
was compared with that of four dasyuromorphian species covering 
a range of craniodental morphologies and feeding ecologies. These 
comprized three extant dasyurids [Dasyurus hallucatus TMM M- 
6921; D. maculatus UNSW Z20; •Sarcophilus harrisii AM10756] and 
one thylacinid {Thylacinus cynocephalus AM1821). Institutional 
abbreviations are Q_MF (Queensland Museum Fossil, Queens- 
land), TMM (Texas Memorial Museum, Austin), UNSW (Uni- 
versity of New South Wales, Sydney) and AM (Australian 
Museum, Sydney). 

We generated 3D finite element models (FEMs) of each skull on 
the basis of computed tomography X-ray (CT) scan data. 
Digimorph (University of Texas; http://www.digimorph.org) was 
the source of CT data of a D. hallucatus skull (0.0784 mm slice 
thickness, 0.0784 mm inter-slice distance). Permission to use 
Digimorph derived CT scan data was granted by Dr. Timothy 
Rowe, Project Director of Digimorph. Other skuUs were scanned 
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in a Toshiba AquUlon 16 scanner (TosliibaMedical Systems 
Corporation, Otawara, Tachigi, Japan) at tlie Mater Hospital, 
Newcastle, NSW (1 mm slice thickness, 0.8 mm inter-slice 
distance, 240 mm field of view). The Australian Museum, 
Queensland Museum and University of New South W ales granted 
the loan of specimens to obtain CT data for this study. 

We used the same specimen CT scans of T. cynocephalus, D. 
maculatus and S. harrisii as examined by Attard et al. [35], but 
constructed the FEMs again using a higher resolution mesh than 
that used by Attard et al. [35] . More specifically, the 3D surface 
meshes which formed the bases of the FEMs were computed using 
the High Quality as opposed to Medium Quality option in Mimics 
(ver. 13.2). The generated surface meshes were then converted to 
FEMs in STRAND7 (ver. 2.4) following previously established 
protocols [31,35,63]. Additionally, 3D objects of T. cynocephalus 
and all dasyurid skulls were exported as separate stl files for the 
creation of interactive 3D pdf documents that can be viewed using 
Adobe Reader (Figure S2, S3, S4, S5). 

Digital Reconstruction of Nimbacinus clicl<soni 

For detailed descriptions of jV dicksoni see Muirhead and Archer 
[14] and Wroe and Musser [18]. Previous analysis has suggested 
that, within the family, the dentition of jV. dicksoni is less derived 
than the recent Thylacinus cynocephalus for at least 12 features, but 
that relatively few cranial specializations in T. cynocephalus 
distinguish the two species. These two taxa share at least three 
cranial features not present in the most generalized thylacinid 
known from significant cranial material, the late Oligocene 
Badjcinus turnhulli [14,18]. 

The skull of JV. dicksoni is well preserved, although some regions 
are absent or damaged. Specifically, some damage/ deformation is 
present at the postorbital processes, frontal, maxillary and nasal 
bones, which are compressed dorsoventraUy (Figure 1). These 
damaged regions were reconstructed according to the morphology 
of surrounding bone regions once the damaged areas had been 
isolated and deleted [64] . Regions of bone that showed only minor 
damage were smoothed to create a coherent surface mesh for later 
solid meshing. 

The right and left dentaries were largely intact but missing the 
superior regions of the coronoid processes, the temporomandib- 
ular joints (TMJ), condyles and angular processes (Figure 1). The 
anterior of the mandible is broken, separating both dentaries. We 
used the right dentary as a basis for reconstruction because its 
dentition was more complete, with only the incisors missing 
(Figure 1). We used a surface mesh of the right dentary of D. 
maculatus to reconstruct posterior regions of the right dentary of jV. 
dicksoni. Dasyums maculatus was chosen as its mandible was most 
similar in shape to that of .A^. dicksoni [1 1], thereby minimizing the 
extent of warping needed (and see below). 

Reconstruction involved scaling the dentary of D. maculatus to 
the same size as that of jV. dicksoni on the basis of skull length 
(condylo-basal). The missing posterior region of the jV. dicksoni 
specimen was then isolated on the D. maculatus specimen and the 
mesh fitted to the existing structure in the mesh oiN. dicksoni using 
Iterative Closest Point (ICP) registration. ICP is an algorithm that 
revises the transformation needed to minimize the distance 
between the points of two partially overlapping meshes. This 
process re-oriented the D. maculatus dentary in accordance with the 
morphology of the jV. dicksoni dentary [65] . The anterior region of 
the D. maculatus dentary was deleted and the posterior region 
'warped' so that overlapping regions of the coronoid process and 
angular process from the D. maculatus mesh matched the existing 
morphology of jV. dicksoni. The manual warping method was used 
as much of the target (fossil) morphology was missing, making it 



A 




Figure 1. Digital reconstruction of Nimbacinus dicltsoni. Original 
(grey) and reconstructed 3D (yellow) in (A) lateral view and (B) dorsal 
view. (C) Pre-processed Finite Element model of N. dicksoni, showing 
jaw musculature represented by trusses. 
doi:10.1371/journal.pone.0093088.g001 

impossible to apply homologous landmarks on both complete [D. 
maculatus) and incomplete fossil (jV. dicksoni) specimens. Procedures 
to warp overlapping skuU regions followed estabhshed protocols 
used by Oldfield et al. and Parr et al. [38,66]. Manual warping 
works by establishing a grid of control points around the complete 
model (note that at this stage the incomplete fossil model has been 
ICP registered with the scaled complete model by matching the 
orientations of the regions of the jaw that are present in both 
models). These control points are then manipulated so that the 
surface morphology of the complete model matches that of the 
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Figure 2. Position of nodes selected on each model to measure von Mises stress. Nodes were selected at equidistant points along the (A) 
mid-sagittal plane, (B) zygomatic arch and (C) mandible to measure the distribution of von Mises stress for each loading case. 
doi:1 0.1 371 /journal.pone.0093088.g002 



target (fossil). This is another variation of Template Mesh 
Deformation [66], but with the template points being the grid 
control points around the complete model rather than homologous 
anatomical points on both models. 

Similarly, the TMJ of the complete {D. maculatus) model was 
warped so that the condyle articulated with and fitted the jV. 
dicksoni cotyle of the cranium, again by using the manual template 
mesh deformation warping method. The left dentary was created 
by mirroring the reconstructed right dentary. These were 
positioned so that the condyles articulated with the cranium, the 
outer surfaces of the lower molars made contact with the inner 
surface of the upper molars, and the tips of the lower canines 
aligned with their 'sockets' in the cranium (see Figure 1). 

It is important to note that the shape of the warp was 
determined by the existing regions of the jV. dicksoni dentary; the 
need for the condyle to articulate with the cotyle to form the TMJ 
and for the coronoid process to fit between the cranium and the 
zygomatic arch. These requirements act as restraints on the warp 
such that the shape of the starting mesh [D. maculatus in this case) is 
not important in the sense that the warping process would always 
end with a similarly shaped posterior region of the mandible 
regardless of which taxon was used. We reiterate that D. maculatus 
was used because it was the most similar in shape [11] and 
therefore required less 'warping'. 

The jV. dicksoni cranium was missing the following teeth: left II- 
4, right II, 3-4, both right and left CI and right LDP2. The 
existing 12 and LDP2 on M. dicksoni were mirrored. All incisors 
were missing from the mandible. Incisors from D. maculatus were 
isolated, scaled and fitted into the empty tooth sockets on jV. 
dicksoni. Figure 1 displays the completed reconstruction of jV. 
dicksoni. 

Finite Element Models 

The assembly of FEMs largely follows previously published 
procedures [30,31,35]. As the skull of jV. dicksoni was not fuUy 
preserved, we were unable to assign multiple material properties to 
the digital reconstruction without introducing additional assump- 
tions. Consequendy, as in most FEA incorporating fossil material 
[30,39,67], all FEMs were homogeneous and assigned a single 
material property for cortical bone (£ = 13.7 GPa, » = 0.3, where £ 
is Young's modulus of elasticity and v is Poisson's ratio) [68] to 
enable direct comparisons between species. Poisson's ratio and 



Young's modulus are fundamental metrics in the comparison of 
stress or deformation for any material when strained elasticaUy, 
including homogeneous materials [69]. Young's modulus is a 
measure of stiffness in the material, whereas Poisson's ratio is the 
negative ratio between transverse strain and longitudinal strain in 
an elastic material subjected to uniaxial stress [70]. 

Each homogeneous model was comprized of four-noded 
tetrahedral elements or 'bricks' (tet4). FEMs for QMF36357, 
AM1821, AM10756, UNSW Z20 and TMM M-6921 were 
comprized of 1564048, 1429714, 1799292, 1402103 and 1956942 
bricks respectively. Tet4 models are theoretically less accurate than 
models comprized of tetlO elements. However, the models used in 
this study are large and any difference in accuracy between results 
from tet4 versus tetlO models will diminish as the number of 
elements is increased. Comparable analyses comparing tet4 and 
tetlO based models much smaller than those used here (<252000 
elements) found differences of <10% [27]. 

Modeling Masticatory Muscle Forces 

Jaw elevators were modeled as seven muscle subdivisions: 
temporalis superficialis, temporalis profundus, masseter superficialis, masseter 
profundus, zjgomatico mandibularis, ptetygoideus intemus and pleiygoideus 
extemus [32]. Proportions used for each jaw muscle division were 
based on muscle mass proportions from a dissected Virginia 
opossum [Didelphis virginiand) [71]. Muscle forces were predicted on 
the basis of maximum cross-sectional areas (CSA) using the 'dry 
skull' method [72]. To improve the accuracy of our CSA 
measurements, we used our FEMs to record the co-ordinates of 
~100 nodes at the perimeter of each muscle cross sectional area 
[73]. The FEM was moved to the correct orientation described by 
Thomason [70] to select nodes outiining the CSA. The node co- 
ordinates were then plotted in plane geometry software, GEUP 5 
(version 5.0.3) and connected to form a multi-sided polygon. The 
area of the polygon was measured to estimate the CSA of each 
major jaw closing muscle. To minimize the incidence of artefacts 
at bite points and muscle origin and insertion areas, surface 
regions at these sites were tessellated using a network of stiff beam 
elements [74]. 

Restraints, Loading Conditions and Scaling 

Dasyurids frequently use a penetrating canine bite to kUl prey 
[44,75-78] which involves the application of a bending load [79]. 
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Table 1. Predicted body mass and masticatory muscle forces for modeled dasyuromorphians. 





Species 


Predicted body mass {kg) Temporalis i 


Tiuscle force (N) Masseteric r 


nuscle force (N) Total muscle force (N) 


Dasyurus hallucatus 


0.78 


67.60 


55.89 


123.49 


Dasyurus macutatus 


2.88 


211.01 


1 78.67 


389.67 


Nimbacinus dic!<soni 


5.25 


282.38 


368.33 


650.71 


Sarcophitus harrisii 


14.20 


300.46 


384.73 


685.19 


Tbylacinus cynocepbalus 


32.49 


706.64 


843.21 


1 ,549.85 



Predicted body mass (l<g) calculated using tlie regression equation for dasyuromorphians provided by Myers [76] based on lower molar row length. Temporalis and 
masseteric muscle forces (N) were calculated based on cross-sectional area [67]. 
doi:l 0.1 371 /journal.pone.0093088.t001 



We simulated bilateral canine biting (intrinsic load) and four 
extrinsic loads to simulate loads generated by struggling prey (axial 
twist, lateral shake, puUback and dorsoventral) for all models using 
protocols described by Attard et al. [35] and following McHenry et 



al. [31]. Extrinsic loads were modeled without applying bite forces 
so as to clearly reveal the different influences of each separate 
loading [31]. A gape angle of 35° was applied in all linear static 
load cases. 



ABC 



F 




123456789 
Location along mandible 



Figure 3. Von Mises stress under a bilateral canine bite in lateral view. The models are subjected to a load applied to both canines, with bite 
force scaled based on theoretical body mass. Species modeled were (A) Dasyurus hallucatus, (B) Dasyurus maculatus, (C) Sarcophilus harrisii, (D) 
Nimbacinus dicksoni and (E) Thylaclnus cynocephalus. White colored regions of the skull represent VM stress above 10 MPa. (F) Distribution of von 
Mises stress was measured from anterior to posterior along the mandible. 
doi:1 0.1 371/journal.pone.0093088.g003 
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VM stress 




Location along mid-sagittal plane 

Figure 4. Von Mises stress under a bilateral canine bite in dorsal view. The models are subjected to a load applied to both canines, with bite 
force scaled based on theoretical body mass. Species modeled were (A) Dasyuws hallucatus, (B) Dasyurus maculatus, (C) Sarcophilus harrisii, (D) 
Nimbacinus dicksoni and (E) Ttiylacinus cynocephalus. White colored regions of the skull represent VM stress above 10 MPa. (F) Distribution of von 
Mises stress was measured from anterior to posterior along the mid-sagittal plane. 
doi:1 0.1 371 /journal.pone.0093088.g004 



A considerable size range exists between specimens considered 
in the present study. The relationship between bite force and body 
mass is negatively allometric [26,80]. To account for differences in 
body mass, a second series of load cases were solved following the 
scaling procedures of McHenry et al. [31]. Here, for each model, 
an estimate of bite force was made based on regression of body 
mass to predicted bite force for dasyuromorphians [z = 0.6998 (log 
jl')+l .8735, where andj = mass (g) and z = bite force at canines (N)] 
[26], with body mass for each specimen predicted using the 
equation based on lower molar row length [log j =— 1.075+ 
3.209(log x), where * = lower molar length (mm), andj)' = mass (g)] 
as presented by Myers [81]. Muscle forces were then scaled for 



each specimen to achieve bite forces predicted on the basis of body 
mass. FEMs were solved using these scaled muscle forces. 
Prediction of bite force based on body mass using the regression 
equation provided in Wroe et al. [26] is close to that which would 
be expected following a 2/3 power relationship, whereby muscle 
force is proportional to area while body mass is proportional to 
volume [82] . The maximum bite force measured in Newtons (N) 
was also estimated for intrinsic loads (Table S2) using FEMs with 
un-scaled, specimen-specific estimated muscle forces (Table S3). 
Three dimensional approaches are likely to be more accurate than 
2D based approaches [83]. 
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Figure 5. Von Mises stress under extrinsic loads in lateral view. The models are subjected to various loads applied to the canines, including a 
(A, E, I, M, Q) lateral shake, (B, F, J, N, R) axial twist, (C, G, K, 0, S) pullback and (D, H, L, P, T) dorsoventral. The force applied was equivalent to 1 00 times 
the animal's estimated body mass for an axial twist, and 10 times the animal's estimated body mass for a lateral shake, pullback and dorsoventral 
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shal<e. Species compared were (A-D) Dasyurus hallucatus, (E-H) Dasyurus maculatus, (l-L) Sarcophilus harrisii, (M-P) Nimbacinus dicksoni and (Q-T) 
nyiacinus cynoceptiaius. White colored regions of the skull represent VM stress above 10 MPa. Distribution of von IVlises (VIVI) stress was measured 
from anterior to posterior along the mandible for a (U) lateral shake, (V) axial twist, (W) pullback and (X) dorsoventral. 
doi:l 0.1 371/journal.pone.0093088.g005 



A H-frame connecting the canines of the upper and lower jaws 
was used to apply extrinsic forces, with forces applied at the center 
of the frame [32,35]. The force (N) applied to extrinsic loads was 
an arbitrary figure, applied for strictly comparative purposes, 
equivalent to 100 times the animal's estimated body mass for an 
axial twist, and 10 times the animal's estimated body mass for a 
lateral shake, puUback and dorsoventral shake [81]. Each 
simulation in which forces are applied with the anterior teeth 
(canines) restrained is a test for the hypothesis that stress wiU be 
highest for species with the longest rostrum. 

Von Mises (VM) stress is a good predictor of failure in ductile 
materials such as bone [84,85] and VM stress is used here as a 
metric for comparison between models following Attard et al. [35]. 
Nodes were selected at equidistant points along the mid-sagittal 
plane, zygomatic arch and mandible (Figure 2) following Attard et 
al. [35] and at each node values were calculated by averaging VM 
stress recorded in the surrounding elements to assess changes in 
stress magnitudes and distributions under different loadings. 

Principal component analysis (PCA) was used to visualize 
differences between species in average VM stress values for 
equidistant nodes along the mid-sagittal plane (N= 10). PCA is an 
ordination technique that summarizes the maximum variation 
among a set of variables on few, uncorrelated axes (principal 
components) [86]. PCA was performed separately on VM stress 
values for each extrinsic load (axial twist, lateral shake, pullback 
and dorsoventral) and for a bilateral canine bite. All VM stress 
values were log transformed prior to PCA. 

Differences in VM stress values between species were compared 
using a Kruskal-Wallis test, which is regarded as a multiple-group 
extension of the Man- Whitney test [87] . Significance values were 
corrected for multiple comparisons using Bonferroni corrections, 
as a conservative approach. 

Results 

The predicted body mass (kg) of each species was generally 
within the expected range for each of the extant species (Table 1). 
Body mass estimates ranged from 0.78 kg for D. hallucatus, up to 
32.49 kg for T. cynocephalus. However, the body mass estimated for 

5. harrisii of 14.20 kg was slightly abo\'e the upper limit observed 
for males (13 kg) [88], possibly because the teeth and skull are 
relatively large in this species. The robust craniodental morphol- 
ogy and relatively large teeth in S. harrisii are probably related to its 
habitual osteophagy, as has been observed in bone-cracking 
carnivorans [89]. To obtain body mass estimates for these taxa 
using simple or multiple regressions adjusted from cranidoental 
variables may lead to overestimates of body mass. Predicted 
maximum muscle forces for jV. dicksoni (65 1 N) were relatively high, 
approaching those of the larger S. harrisii (685 N) (Table 1). 

Thylacinus cynocephalus displayed comparatively high levels of VM 
stress in the cranium and mandible for most simulations (Figure 3- 

6, S3). This is consistent with results of Attard et.al. [35]. Dasyurus 
hallucatus showed relatively high levels of stress in the posterior of 
the mandible for a canine bite (Figure 3A), and along the ventral 
surface of the ramus for most extrinsic loads (Figure 5A-C). 

The regions of highest stress along the dentary of jV! dicksoni were 
located at the coronoid fossa and condylar process (Figure 3D). 
These regions of peak stress may be in part artifacts of 
reconstruction. Otherwise the dentary of N. dicksoni revealed 



similar stress patterns for a bilateral bite to D. maculcUus (Figure 3). 
The distribution of stress for M. dicksoni in the cranium in response 
to a bilateral bite was intermediate between S. harrisii and D. 
maculatus (Figure 4). The magnitudes of stress along the mid- 
sagittal plane of jV. dicksoni were slightly higher than for S. harrisii 
and lower than for D. maculatus (Figure 4F). 

The highest stress in the cranium occurred at the zygomatic 
arch for all species in response to a bilateral canine bite (Figure 4). 
Von Mises stress along the zygomatic arch during a bilateral 
canine bite gradually increased posteriorly in S. harrisii, while stress 
peaked at node 3 in T. cynocephalus followed by a gradual decrease 
posteriorly (Figure S6A). The three other species displayed two 
peaks in stress along the zygomatic arch during a bilateral canine 
bite; one at the middle and the other at the posterior region of the 
zygomatic arch. Thylacinus cynocephalus was the only species to show 
two distinct peaks in stress for a bilateral bite along the mid-sagittal 
crest (Figure 4F). These stress points occurred at the temporal 
ridge and at the most narrowed region of the nasal (Figure 4E). 
Von Mises stress measured along the mid-sagittal crest for a 
bilateral bite revealed one point of peak stress halfway along the 
frontal of Z). maculatus, S. harrisii and jV! dicksoni and at the temporal 
ridge for D. hallucatus (Figure 4). 

Stress was quite evenly distributed along the dentaries for all 
species in response to lateral shaking and axial twisting, with the 
exception of D. hallucatus, wherein stresses peaked anteriorly 
(Figure 5U-V) resulting in significandy higher VM stress values for 
that species compared to all others (x^ = 32.87, P<0.0001). An 
axial twist resulted in much higher lc\ (;ls of str(\ss along the mid- 
sagittal crest for D. hallucatus compared to all other species, and 
peaked at the anterior of the nasal and at the frontal (Figure 6B). 
Sarcophilus harrisii and T. cynocephalus showed higher levels of stress 
along the mandible for a pullback and dorsoventral shake than 
other species included in this study (Figure 5W— X). Comparisons 
of mandible VM stress values revealed significant differences 
between species after Bonferroni correction for both pullback 
(X^ = 33.28, P<0.001) and dorsoventral shake (x^ = 35.61, P< 
0.0001), with the exception of T. cynocephalus and S. harrisii. Two 
points of peak stress were apparent along the dentary for T. 
cynocephalus in these two simulations; one at the most anterior 
point, and the second at the coronoid fossa. Stress distribution 
along the dentary of S. harrisii followed a similar trend for a 
pullback and dorsoventral shake; peaking at the ramus inferior to 
Ml then gradually decreasing posteriorly. 

PCA results for mid-sagittal node VM stress values (Figure 7) 
showed that a high proportion of variance could be explained in 
all cases by two Principal Component (PC) axes (>85%). These 
plots provide an appreciation of interspecific difiFerences across all 
10 mid-sagittal nodes and bite simulations. PCA results indicate 
that the main axes of interspecific variance for all bites were 
explained by either nodes 1 and/or 7-10. Thylacinus cynocephalus 
and S. harrisii difiFered significantly for VM stress values under a 
bilateral bite at the canines (X^ = 12.95, P= 0.04) and PCA results 
indicated separation of those two species along PCI (60.6%), 
which largely explained variance in node 8 and node 10 
(Figure 7A). PCI for a bilateral canine bite revealed close 
similarities between jV. dicksoni and D. maculatus, whereas PC2 
(28.9%) clearly separated jV. dicksoni from D. maculatus and reflected 
differences in node 7 (as seen in Figure 4). 



PLOS ONE I www.plosone.org 



8 



April 2014 | Volume 9 | Issue 4 | e93088 



Prey Size Preference in Nimbacinus dicksoni 



A B C D 




Q R S 




123456789 10 123456789 10 

Location along mid-sagittal plane Location along mid-sagittal plane 



Figure 6. Von Mises stress under extrinsic loads in dorsal view. The models are subjected to various loads applied to the canines, including a 
(A, E, I, M, Q) lateral shake, (B, F, J, N, R) axial twist, (C, G, K, 0, S) pullback and (D, H, L, P, T) dorsoventral. The force applied was equivalent to 1 00 times 
the animal's estimated body mass for an axial twist, and 10 times the animal's estimated body mass for a lateral shake, pullback and dorsoventral 
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shake. Species compared were (A-D) Dasyurus hallucatus, (E-H) Dasyurus maculatus, (l-L) Sarcophilus harrisii, (M-P) Nimbacinus dicl<soni and (Q-T) 
Thylacmus cynocephalus. White colored regions of the sl<ull represent VM stress above 10 MPa. Distribution of von Mises (VM) stress was measured 
from anterior to posterior along the mid-sagittal plane for a (U) lateral shake, (V) axial twist, (W) pullback and (X) dorsoventral. 
doi:1 0.1 371 /journal.pone.0093088.g006 



For a lateral shake, PCI (51.0%) explained change at nodes 1 
and 10 and separated T. cynocephalus from jV. dicksoni (Figure 7B), as 
also seen in Figure 6U. PC2 (42.0%) for a lateral shake revealed 
difierences among all species for nodes 7 and 10, with T. 
cynocephalus closely resembling jV. dicksoni compared to other species 
(Figure 7B). Interspecific differences were not significant for a 
lateral shake after Bonferroni correction, and before correction 
those distinguished T. cynocephalus from JV. dicksoni, D. maculatus and 
S. harrisii (x-' = 9.27, /"= 0.01-0.03). PCA for an axial twist 



(Figure 7C) showed that the extremes of PC 1 (66. 1 %) were 
delimited by S. harrisii and D. hallucatus, and the remaining taxa 
were located in between those two species. PC 1 mainly accounted 
for difiFerences among species for node 1 values (as seen also in 
Figure 6V), which were high for D. hallucatus and low for S. harrisii, 
whereas PC2 (26.1%) largely summarized node 10 values. Pairwise 
comparisons of VM stress values were not significant for an axial 
twist, however before Bonferroni correction, difiFerences between 
S. harrisii and D. hallucatus (% =10.73, P= 0.03), jV. dicksoni 
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Figure 7. Principal components one and two of von Mises stress along mid-sagittal plane. Results of principal components analysis to 
compare stress among species along mid-sagittal nodes for each loading case, including (A) bilateral canine bite, (B) lateral shake, (C) axial twist, (D) 
pullback and (E) dorsoventral. Key to symbols: (pink square) Thylacmus cynocephalus; (red triangle) Nimbicinus diclisoni; (orange diamond) Sarcophilus 
harrisii; (green circle) Dasyurus maculatus; (purple triangle) Dasyurus hallucatus. 
doi:10.1371/journal.pone.0093088.g007 
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(X^= 10.73, P=0.01) and T. cymcephahs (x^= 10.73, P=0.04) 
could be distinguished from one another. For a pullback bite, PC 1 
(50. 1 %) mainly explained interspecific differences for node 10, and 
PC 2 (43.ry()) explained change at nodes 1 and 8 (Figure 7D). VM 
stress values for T. cjnocephalus were different from those for D. 
hallucatus (x^= 14.37, P=0.02) and D. maculatus (x^= 14.37, 
P=0.04) for a pullback bite. For the dorsoventral shake, PCI 
explained 60.2% of variance and reflected differences between T. 
cynocephalus and S. harrisii, whereas values across nodes were more 
similar among the remaining three species, located toward the 
middle of PCI (Figure 7E). PC2 (31.4%) separated D. hallucatus 
from T. cynocephalus, and pairwise comparisons revealed node 
values to be different between those two species {y^ = 14.78, 
P=0.01). 

Discussion 

Differences in biomechanical performance between the three 
extant dasyurids included in this study appear consistent with their 
respective known feeding behaviors. Dasyurus hallucatus showed 
comparatively higher levels of stress in most simulations than S. 
harrisii and D. maculatus. Dasyurus hallucatus eats invertebrates and 
other relatively small prey [60-62], which may not require 
adaptation to sustain the fuU range of extrinsic loads simulated 
here. This species shows particularly high VM stress in axial 
twisting, especially in contrast to S. harrisii However, it performs 
relatively well under pull-back loading, which may be linked to a 
capacit)' for pulling invertebrates from the ground. Observational 
studies on wild D. hallucatus will be required to confirm the 
functional role of their skull in prey acquisition. Future work on 
the comparative musculoskeletal anatomy and collection of in vivo 
or ex vivo biomechanical data of the extant species would likely 
improve the predictive power of current bite force and muscle 
force estimations. Overall consistencies found between known prey 
size and biomechanical performance for extant dasyuromorphians 
underscore the potential value of projections based on compar- 
ative FEA for extinct/fossil taxa. 

Our comparative biomechanical modeling of dasyuromorphian 
skulls suggests considerable differences in predatory behaviors 
between the two thylacinids considered here. Our 3D based results 
indicate that the Oligocene to Miocene jV. dicksoni had a high bite 
force for its size, comparable to that of extant dasyurids known to 
take relatively large prey, D. maculatus and S. harrisii [48,60]. In 
light of similar levels of 'carnassialization' (development of 
relatively long, high amplitude vertical shearing crests) in the 
cheektooth dentition with D. maculatus, and a lack of obvious dental 
specialization consistent with regular bone-cracking, our results 
suggest a predominandy carnivorous diet for N. dicksoni that may 
have included relatively large prey. Dasyurus maculatus are 
opportunistic hunters, varying their diet in response to environ- 
mental disturbances and short-term fluctuations in prey abun- 
dance [54,56]. They will prey on vertebrate species up to and 
sometimes exceeding their own body mass. Prey includes 
bandicoots, smaller dasyurids, possums, smaller macropodoids, 
snakes, lizards, birds and frogs, as well as invertebrates. Potential 
prey for a fox-sized thylacinid living in the closed forest 
communities of Riversleigh likely included many small to 
medium-size birds, frogs, lizards and snakes, as well as a wide 
range of marsupials, including bandicoots (peramelemorphians), 
dasyurids (dasyuromorphians), kangaroos (macropodoids), thingo- 
dontans (yaJkiparidontians), marsupial moles (notoryctemorphians) 
and wombats (vombatoids) [12]. 

Although our FEA results for N. dicksoni suggest a capacity to kill 
prey approaching or exceeding its own body mass, its prey range 



may have been limited by competition with sympatric carnivores. 
The extent of niche overlap and competition within this ancient, 
medium-large sized carnivore community may have been partially 
alleviated by occupying different habitats and specializing in 
different hunting strategies. The recovery of a near complete 
skeleton of jV. dicksoni [14] will provide further information on the 
locomotion and predatory behavior based on postcranial material; 
for example, was N. dicksoni as arboreal as the extant D. macuhtusi 

Diffcrc-nccs in mechanical performance suggest that T. 
cynocephalus is unusual relative to other dasyuromorphians, 
including, jV. dicksoni, as indicated by distinctly higher VM stresses 
than all other species in response to each loading case. Thylacinus 
cynocephalus, in contrast to jV. dicksoni, has completely lost the 
metaconid on the lower molars and has a proportionately much 
larger postmetacrista on the upper molars. On the basis of 
traditional beam theory we predicted that taxa with longer rostra 
would exhibit higher stress [90] , as evident in the long-snouted T. 
cynocephalus relative to shorter-snouted dasyuromorphians. Differ- 
ences between T. cynocephalus and other species were also 
significant for three out of five simulations examined after 
conser\'ative Bonferroni correction for multiple testing. These 
results further support the contention by Attard et al. [35] that 
niche breadth in T. cynocephalus may have been more limited and 
that it likely preyed on relatively small to medium-sized vertebrates 
such as wallabies, possums and bandicoots. 

Although measures of skull performance in response to forces 
imposed by struggling prey revealed closer similarity between the 
fossil thylacinid N. dicksoni and large extant carnivorous dasyurids, 
than with T. cynocephalus, there were differences. Our reconstruc- 
tion suggests that the TMJ was more elevated in N. dicksoni than in 
D. maculatus, and higher relative to the height of the cheektooth 
row. The TMJ is a complex joint and is important for occlusion 
and mastication [91,92]. The position of the TMJ can inffuence 
bite strength and muscle activation [93] . The position of the TMJ 
Eilong the anterior-posterior axis tends to lie closer to the plane of 
the tooth row in carnivorous taxa [1 1]. Conclusive determination 
of the precise position and morphology of the TMJ in N. dicksoni 
must await the discovery of more complete cranial material. 

Morphological evidence from past studies further demonstrates 
diversit)' within this family. The smallest thylacinid, Muribacinus 
gadiyuli, is believed to have fed on relatively small vertebrates and 
invertebrates because it lacks a number of dental features present 
in large prey specialists (e.g., robust protoconids and brachyce- 
phalization) such as similarly sized D. maculatus [17]. The variety of 
feeding behaviors among thylacinids may have helped facilitate 
their co-existence within different ecological niches that were later 
filled by diversifying carnivorous dasyurids. 

Supporting information 

Figure SI Phylogenetic tree of dasyuromophians inves- 
tigated in this study. One of several recent assessments of the 
phylogenetic relationships of dasyuromorphians, including taxa 
that have been examined in this study (Wroe & Musser 2001). 
(TIF) 

Figure S2 Interactive 3D pdf showing the digitally 
segmented cranium and mandible of TTiylacinus cyno- 
cephalus. 

(PDF) 

Figure S3 Interactive 3D pdf showing the digitally 
segmented cranium and mandible of Sarcophilus harri- 
sii. 

(PDF) 



PLOS ONE I www.plosone.org 



11 



April 2014 | Volume 9 | Issue 4 | e93088 



Prey Size Preference in Nimbacinus dicksoni 



Figure S4 Interactive 3D pdf showing the digitally 
segmented cranium and mandible ot Dasyurus macula- 
tus. 

(PDF) 

Figure S5 Interactive 3D pdf showing the digitally 
segmented cranium and mandible o£ Dasyurus halluca- 
tus. 

(PDF) 

Figure S6 Von Mises stress along zygomatic arch for all 
loading cases. Distribution of von Mises (VM) stress was 
measured from anterior to posterior along the zygomatic arch for 
a (A) bilateral canine bite, (B) lateral shake, (C) axial twist, (D) 
puUback and (E) dorsoventral. 
(TIF) 

Table SI Temporal and geographic distribution of 
thylacinid species. Abbreviations: Aust, Australian mainland; 
E., Early; L., Late; M., Middle; Mio, Miocene; NG, New Guinea; 
NT, Northern Territory; Oligo, Oligocene; Plio, Pliocene; Qld, 
Queensland; Tas, Tasmania. 
(PDF) 

Table S2 Maximum bite forces (N) for un-scaled 
homogeneous models for a bilateral canine bite. 

(PDF) 
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